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1. Introduction 

Many experimental data analysis and model validation consist in solving a so-called 
inverse problem [1]. Typically, one wishes to determine the distribution of a hidden 
variable knowing the distributions of measured variables and the transform function that 
connects both types of variable^ As a distribution is, by definition, a positive function, 
the algorithm used to solve this kind of inverse problem must converge towards the 
optimum positive solution. A very popular algorithm used to solve this problem is the 
Hanson and Lawson NNLS (Non- Negative Least Squares [2j§|. However, this algorithm 
consists in minimizing the residue between the measured observable distribution and 
the theoretical distribution (i.e. the hidden distribution times the transform function). 
In most actual situations, the measured distributions and the transform matrices are 
affected by fluctuations and uncertainties. In such a case, what has to be minimized 
is not a residue but a chi-square. In the following, we show how these different kinds 
of uncertainties may be taken into account in the minimization. We introduce a new 
algorithm NNLC (Non-Negative Least Chi-square) based on NNLS which allows to 
handle uncertainties. In the second part, we apply both algorithms to a practical 
example. It will be shown how the new protocol improves the gamma hit location 
inside germanium detectors. 

2. Solving of the inverse problem 
2.1. Formalism 

We consider a discreet inverse problem: 

Ax = b, (1) 

where A is the transform matrix, x is the vector of unknowns and b is the (measured) 
vector of observations. When x represents the unknown distribution of a variable then 
each of its components has to be non- negative [3] . On the other hand, the components of 
the transform matrix and of the observation matrix may be negative (this is the case if 
these distributions are detector signals). Most of the time, in real- life problems, no exact 
solution exists and what is searched is the non-negative solution x*^"^ that minimizes the 
residue: 

= ^(6i-6f )2 = ||b- Ax^°'f (2) 

i 

with b^°i = Ax^°^ (3) 

In the following the j**^ columns of matrix A will be considered as a vector and noted 
a,. 

I For example, the determination of energy spectrum of the particles entering a detector from the shape 
of the delivered signals. In this case, the transform function is the response function of the detector. 
§ NNLS is implemented in the Matlab environment as Isqnonneg. 
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2.2. The NNLS algorithm 

Such linear systems can be solved using the NNLS method. This fast iterative algorithm 
converges towards a positive x*^"' vector with a maximum of null components. Indeed, 
the initial guess for the solution is the null-vector. The main loop of the algorithm 
consists in adding one minimization component to the system solving until the stopping 
condition is reached. Thus, at the k^^ iteration, the size of the matrix to invert is k x k. 
At each iteration, the new component added to the system is the one that maximizes the 
gradient of the residue. The stopping condition is reached when the absolute value of this 
gradient is lower than a given tolerance. When the inversion gives negative components, 
an inner loop adds a kernel vector to the solution until it becomes non- negative. 

3. Chi-square minimization 

3.1. Formalism 

The mere least square optimum solution of Eq. ([T| gives often an unsatisfactory solution 
since the components of the right hand side vector and of the transform matrix may be 
affected by fluctuations (noise, statistical fluctuations, . . . ), uncertainties or biases. In 
the most usual cases, what has to be minimized is not the residue, but a chi-square: 



where cxb and cryysoi are the standard deviations of the fluctuations on b and b'^" . The 
first term is usually connected to the limitations of the detection system, the second 
term results from the uncertainties on the transform matrix, Eq. ([s]). 

We first consider the case when the transform matrix is not affected by uncertainties 
(o"]-,soi = 0). Normalizing b and A by the standard deviation of the fiuctuations of the 
observations: b[ = hi/ab^ and a^^ = aij/abi, one obtains = X]j '^ij ^T^- Hence, the 
chi-square can be written residue: 



Thus, least square minimization algorithms can also be used to solve this maximum 
likelihood problem. 

However, in many cases, the uncertainties (characterized by the standard deviations 
(Ja^.) on the transform matrix components cannot be neglected. These uncertainties 
may result from the fact that the transform matrix was obtained from a measurement 
or through a theoretical model that entails simplifications or hypotheses. Considering 
Eq. (|3]), the uncertainties on the transform matrix induce a second denominator term 
in the chi-square: 




(4) 




(5) 




(6) 
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Therefore, when the transform matrix is uncertain, the chi-square cannot be turned 
directly into a residue as in Eq. ([5|, since cXbsoi depends on x*^"' which is the unknown to 
be determined. However, if the optimization algorithm used to solve the linear system 
works iteratively as it is most often the case (NNLS is iterative), the x*^°^ obtained at 
iteration k — 1 can be used to calculate the fluctuations at iteration k. Hence, at each 
iteration we apply the following normalizations: 



so that: 



I 



- -k.^Y.<A'"^ ' (7) 



(fc) Oi^ 
"'ij (k) ' 

cr- 



IbW-A^'^^xWf . (10) 
The chi square is finally transformed into a residue. 



3.2. The NNLC algorithm 

In order to generalize the NNLS algorithm to inverse problems affected by uncertainties, 
we introduce the NNLC algorithm (Non-Negative Least Chi-square). As we have 
seen writing Eq (|5]), when only the right hand side observation vector is affected by 
uncertainties, then the normalization of b and A can be made prior to the iterations. 
In the general case, the normalizations of Eqs. ([^|9| have to be included in the main 
loop. Some caution has to be taken when the uncertainties on A are large with respect 
to the uncertainties on b. In this case, the normalization factor a^'^^ may vary a lot 
from one iteration to the next one, Eq. ([t]). This means that the linear system changes 
at each iteration and that the minimization may not converge, or converge very slowly. 
When this effect appears, it can be solved by allowing a maximum variation of the 
standard-deviations of 10% at each step. 

In the following, we apply this new algorithm to the problem of locating the hits 
of gamma-rays inside a detector. 



4. Application to the location of the interactions of a gamma-ray into a 
germanium crystal 

4.1. Signal decomposition as an inverse problem 

A gamma-ray which enters a germanium detector may interact one or more times with 
the electrons of the crystal. In order to use this property for gamma detection, an 
electric field is applied to the crystal, so that each gamma-electron interaction (hit) 
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provokes electron and hole cascades towards the anode and the cathode. The motions 
of the charges induce charge-signals at the electrodes. The shapes of the signals depend 
on the locations of the hits. When several hits occur simultaneously, the resulting signal 
is the sum of the signals induced by the individual hits. Moreover, the amplitude of the 
signal is proportional to the energy deposited by the gamma. Finding the locations and 
the energies of the hits allows to measure the energy and the direction of the gamma-ray 
[3]. This is a widely used technique for modern gamma detection [H El El IE] ■ 

The problem of determining the locations of the hits from the shape of the resulting 
signal [9] is a typical inverse problem. The observable vector is the detected signal, the 
vector of unknowns contains the energy deposit on each point of the crystal and the 
transform matrix gives the shape of the signals for every point in the crystal. Of course, 
in order to keep finite size A and x matrices, the volume of the crystal is discretized 
into a finite number of voxels (in our case, we consider 2x2x2 mm voxels for a crystal 
volume of about 400 cm^). Each column of the transform matrix is then the signal 
corresponding to a hit deposing a unit energy in the j*^ voxel, and xj is the energy 
actually deposited. Due to the additive property of the signals, signal decomposition 
corresponds to the solving of Eq. ([T]). 




Figure 1. AGATA germanium crystal and its capsule. The hole along the axis is the 
anode. The outer surface is covered by 36 cathodes. 

The AGATA [TO] germanium crystal used in this analysis is shown in Fig. [l} The 
single anode is along its central axis. The outer surface is covered by six slices of six 
cathodes, thus the volume of the crystal is divided into 36 electrical segments. When a 
hit occurs in a given segment, the corresponding cathodes measures a net charge signal 
and the neighboring segments measure transient signals. All these signals, concatenated 
the one after the other, are used to form the b and the vectors (see Fig. [2]). 

The SLj signals forming the transform matrix are obtain using the MGS simulation 
code [IT] (they may also be obtained using a scanning device [I2l 113]). 

4-2. Uncertainties 

The main fluctuations on the detector signals are due to the electronic noise. This 
noise does not depend on the signal sample and corresponds to a standard-deviation 
= CTnoise ~ 3 keV. Moreovcr, the continuous signals delivered by the detector are 
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Figure 2. Example of signal concatenation used to build the b and the vectors. 
The 9 signals belong to the hit segment (fourth signal) and to its neighbors. 



time-discretized. As the discretizer is triggered when the noisy signal crosses a given 
threshold, the resulting discretized signals can be slightly time shifted p3l [T5] (in our 
case, all the segment signals from a given event are translated by the same time shift). 
Typically, for tg = 10 ns samples, the shift At is of the order of some nanoseconds. The 
influence of the time shift on the amplitude in a given sample is illustrated in Fig. |3j 



1.5- 




i-l i i+l 
time index 

Figure 3. The time shift At induces an amplitude shift Ab proportional to the slope 
s. 



For small time jitters, the effect on the amplitude is proportional to the signal slope 
and to the time shift: A6j = — Sj At. The slope of the signal can be estimated by: 

bi+i — hi_i 

= —itr~ ' ^''^ 

where ts is the sampling duration. 
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The slope calculated using the previous equation for the detected signals would be 
excessively affected by the noise (the fluctuations increase with the degree of derivation). 
Thus it is necessary to consider that the detected signal is the time reference, and that 
the transform matrix signals are time shifted with respect to the detected signal. This 
way, the slopes are calculated from smooth simulated signals. The fluctuations on the 
transform matrix elements are: 



— a,- 



2t 



O'At 



(12) 



where a/^t is the standard deviation of the time shifts [HI 

Finally, from Eq. (|4]), the denominator of the chi-square reads: 
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This expression of the standard deviation is declared to the NNLC algorithm. 
4-3. Results 

In order to compare the NNLS and NNLC algorithm results in this typical situation 
where both the observation vector and the transform matrix are altered by uncertainties, 
we simulate a large number of hits, at known random positions inside the crystal, 
for different values of the noise and of the time jitter standard deviation. We then 
apply both algorithms to calculate x*^"^ and compare the precisions on the resulting hit 
locations. 




Figure 4. Mean error on the location of the hits (averaged over 350 hits per point). 
The open dot curves are obtained using NNLS, the black dot curves are obtained using 
NNLC. The dotted lines correspond to 30 keV energy deposits, the thin lines to 300 
keV deposits and the bold lines to 3 MeV deposits. The noise on the signals is 3 keV 
and the sample time is 10 ns. 



The precision of the location of the hit is presented in Fig. |4] as a function of the 
time jitter. The lower curves correspond to 3 MeV energy deposits, the middle curves 
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to 300 keV deposits and the upper curves to 30 keV energy deposits. The error is lower 
for higher energies since the signal-to- noise ratio is larger. When there is no time jitter, 
the fluctuations are due to the noise, thus their standard deviation does not depend 
on the sample. In this case, the denominator can be factorized from the chi-square: 
= (1/^noise) ~ ^i"')^ ^ ^'i^- ( 4p3 ) and the chi-square and the residue 

minimizations give the same result. This is no longer true when the signals are affected 
by a time jitter even for small shifts (much lower than 10 ns, the sample duration). 
When the effect of the noise becomes negligible with respect to the effect of the time 
shift, the error on the position is almost independent on the energy deposit. For positive 
time jitters, the NNLC algorithm always gives better results than the NNLS algorithm 
and, as expected, the difference increases with the time shift. 



5. Conclusions 



This paper introduces a new method to account for the fluctuations and the uncertainties 

that may affect the transform function of an inverse problem as well as the measured 

observation vector. This method is embodied in the NNLC algorithm which is an 

evolution from the NNLS algorithm that replaces residue minimization by chi-square 

minimization. This new method was applied to the location of gamma-ray interactions 

inside germanium detectors. It was shown that, whatever the signal-to- noise ratio and 

the time jitter, the precision on the location is always better using NNLC. We believe 

that this very general method has a large number of actual applications. 
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